Lecture 9

Spatio-temporal models

Sara Martino
sara.martino@ntnu.no

Dept. of Mathematical Science, NTNU

Janine Illian
Janine.Illian@glasgow.ac.uk

University of Glasgow

Jafet Belmont
Jafet.Belmont@glasgow.ac.uk

University of Glasgow

Motivation

  • Many real-world data sets vary in both space and time.

    • Examples: rainfall, temperature, air pollution, crop yield, disease spread.

  • We often want to model how similarity (correlation) changes across both dimensions.

The good news 👍

  • the framework works with spatio-temporal data as well

The bad news 👎

  • fully spatio-temporal models are complex

  • model fitting can take a long time

  • different types of spatio-temporal data structures lead to different types of complexities

Complexity of spatio-temporal models

  • additional dependencies: in space and time

  • spatial and temporal behaviour can be independent on each other, or dependent: i.e. properties of the spatial structure vary through time (or vice versa)

Several options are available in inlabru

  • Areal Model

  • Geostatistical and Point Process models

The Spatio-Temporal Framework

We define a stochastic process:

\[ Z(s,t), \quad s \in \mathcal{S} \subset \mathbb{R}^d, \quad t \in \mathcal{T} \]

The covariance function:

\[ C((s_1,t_1),(s_2,t_2)) = \text{Cov}[Z(s_1,t_1), Z(s_2,t_2)] \]

Our goal: specify \(C(\cdot)\) to capture spatial, temporal, and joint dependencies.

Separable and non-separable models

Separable models \[ C((s_1,t_1),(s_2,t_2)) = C_S(s_1,s_2) \times C_T(t_1,t_2) \]

Interpretation:

  • Space and time act independently.

  • Spatial correlation not dependent on time lag.

Advantages: 😀

  • Simpler estimation.

  • Fewer parameters \(\rightarrow\) Easier computation.

Disadvantages: 😔

  • Unrealistic for evolving or propagating phenomena.

Non-Separable models

Separable and non-separable models

Separable models \[ C((s_1,t_1),(s_2,t_2)) = C_S(s_1,s_2) \times C_T(t_1,t_2) \]

Interpretation:

  • Space and time act independently.

  • Spatial correlation not dependent on time lag.

Advantages: 😀

  • Simpler estimation.

  • Fewer parameters \(\rightarrow\) Easier computation.

Disadvantages: 😔

  • Unrealistic for evolving or propagating phenomena.

Non-Separable models

\[ C((s_1,t_1),(s_2,t_2)) = f(|s_2-s_1|,|t_2-t_1|) \] Advantages: 😀

  • More realistic for dynamic physical processes.

  • Can model propagation or decay that varies over time.

Disadvantages: 😔

  • More parameters → complex estimation.

  • Need more data for estimation

  • Interpretation may be less straightforward.

We here focus on separable models

Areal Data

Spatio temporal models for areal data

\(\frac{\text{Number of cases } Y_{st}}{\text{Expected numer of cases }E_{st}}\) in Ohio from 1968 to 1988 (simulated data!)

Space

Spatio temporal models for areal data

\(\frac{\text{Number of cases } Y_{st}}{\text{Expected numer of cases }E_{st}}\) in Ohio from 1968 to 1988

Time

Spatio temporal models for areal data

\(\frac{\text{Number of cases } Y_{st}}{\text{Expected numer of cases }E_{st}}\) in Ohio from 1968 to 1988

Time

Spatio temporal models for areal data

The observation model

\[ Y_{st}|\lambda_{st}\sim\text{Poisson}(E_{st}\lambda_{st}) \]

We are going to see 3 different models for the linear predictor \(\eta_{st} = \log \lambda_{st}\)

  1. Model 1: \(\eta_{st} = \beta_0 + u_s + v_t\) with \(u_s\sim\text{CAR}(\tau)\)

    Here the time effect \(v_t\) is constant w.r.t. space and the space effect \(u_t\) is constant w.r.t. time.

Model 1 - Implementation

cmp_1 = ~ Intercept(1) +
  space(id_area, model = "besag", graph = g, scale.model = T) +
  time(year, model  = "rw2", scale.model = T)

lik = bru_obs(formula = Y~.,
              data = out,
              family = "poisson",
              E = E )
fit_1 = bru(cmp_1, lik)

Model 1 - Implementation

cmp_1 = ~ Intercept(1) +
  space(id_area, model = "besag", graph = g, scale.model = T) +
  time(year, model  = "rw2", scale.model = T)

lik = bru_obs(formula = Y~.,
              data = out,
              family = "poisson",
              E = E )
fit_1 = bru(cmp_1, lik)

Spatio temporal models for areal data

The observation model \[ Y_{st}|\lambda_{st}\sim\text{Poisson}(E_{st}\lambda_{st}) \]

We are going to see 3 different models for the linear predictor \(\eta_{st} = \log \lambda_{st}\)

  1. Model 1: \(\eta_{st} = \beta_0 + u_s + v_t\) with \(u_s\sim\text{CAR}(\tau)\)

    Here the time effect \(v_t\) is constant w.r.t. space and the space effect \(u_t\) is constant w.r.t. time.

  2. Model 2: \(\eta_{st} = \beta_0 + u_{st} + v_t\) with \(u_{st}\sim\text{CAR}(\tau),\ i = 1,\dots,T\)

    Here the spatial fields are different for each year, they are independent on each other.

Model 2 - Implementation

cmp_2 = ~ Intercept(1) + 
  space(id_area, model = "besag", graph = g, scale.model = T,
                               replicate = id_time) +
  time(year, model  = "rw2", scale.model = T)

fit_2 = bru(cmp_2, lik)

Model 2 - Implementation

cmp_2 = ~ Intercept(1) + 
  space(id_area, model = "besag", graph = g, scale.model = T,
                               replicate = id_time) +
  time(year, model  = "rw2", scale.model = T)

fit_2 = bru(cmp_2, lik)

Spatio temporal models for areal data

The observation model

\[ Y_{st}|\lambda_{st}\sim\text{Poisson}(E_{st}\lambda_{st}) \]

We are going to see 3 different models for the linear predictor \(\eta_{st} = \log \lambda_{st}\)

  1. Model 1: \(\eta_{st} = \beta_0 + u_s + v_t\) with \(u_s\sim\text{CAR}(\tau)\)

    Here the time effect \(v_t\) is constant w.r.t. space and the space effect \(u_t\) is constant w.r.t. time.

  2. Model 2: \(\eta_{st} = \beta_0 + u_{st} + v_t\) with \(u_{st}\sim\text{CAR}(\tau),\ i = 1,\dots,T\)

    Here the spatial fields are different for each year, they are independent on each other.

  3. Model 3: \(\eta_{st} = \beta_0 + u_{st} + v_t\) with \[ u_{st} = \phi u_{st-1} + w_{st}, \text{ with } w_{st}\sim\text{CAR}(\tau),\ i = 1,\dots,T \]

Here we are adding more temporal dependence in the spatial term, using a AR1 type model Here the spatial fields are different for each year, they are independent on each other.

Model 3 - Implementation

cmp_3 = ~ Intercept(1) + 
  space(id_area, model = "besag", graph = g, scale.model = T,
                               group = id_time, control.group = list(model = "ar1")) +
  time(year, model  = "rw2", scale.model = T)

fit_3 = bru(cmp_3, lik)

Model 3 - Implementation

cmp_3 = ~ Intercept(1) + 
  space(id_area, model = "besag", graph = g, scale.model = T,
                               group = id_time, control.group = list(model = "ar1")) +
  time(year, model  = "rw2", scale.model = T)

fit_3 = bru(cmp_3, lik)

Spatio temporal models for areal data

The observation model

\[ Y_{st}|\lambda_{st}\sim\text{Poisson}(E_{st}\lambda_{st}) \]

  1. Model 1: \(\eta_{st} = \beta_0 + u_s + v_t\) with \(u_s\sim\text{CAR}(\tau)\)

  2. Model 2: \(\eta_{st} = \beta_0 + u_{st} + v_t\) with \(u_{st}\sim\text{CAR}(\tau),\ i = 1,\dots,T\)

  3. Model 3: \(\eta_{st} = \beta_0 + u_{st} + v_t\) with \[ u_{st} = \phi u_{st-1} + w_{st}, \text{ with } w_{st}\sim\text{CAR}(\tau),\ i = 1,\dots,T \]

Compare

  • Differences between Model 2 and 3 are not evident when estimating the past, but they would give different predictions

  • Model 3 in this case seem to better fit the data (more about these measures later)

    model      DIC     WAIC
1 Model 1 14013.57 17381.01
2 Model 2 12346.99 12240.29
3 Model 3 11943.31 11923.44
                     mean    sd 0.025quant 0.975quant
Precision for space 13.05  1.28      10.68      15.71
GroupRho for space   0.89  0.01       0.86       0.91
Precision for time  64.79 27.36      26.11     131.75

Spatio temporal models for areal data in inlabru

  • Models can be stuck together using either

    • replicate - the different slices are independent but share the hyperparameter
    • group - different dependence structures are implemented

Types of group model

names(inla.models()$group)
[1] "exchangeable"    "exchangeablepos" "ar1"             "ar"             
[5] "rw1"             "rw2"             "besag"           "iid"


Note The group and replicate features can be used for more than space-time modeling!

Continuous space models

Continuous space models

(Simulated data)

Spatio temporal models for areal data

The observation model

\[ Y_t(s)|\mu_t(s)\sim\mathcal{N}(\mu_t(s),\sigma_y^2); \qquad \eta_t(s) = \mu_t(s) \]

We are going to see 3 different models for the linear predictor \(\eta_{st}\)

  1. Model 1: \(\eta_t(s) = \beta_0 + u(s) + v_t\) with \(u(s)\sim\text{GF}(\rho,\sigma_u)\)

    Here the time effect \(v_t\) is constant w.r.t. space and the space effect \(u(s)\) is constant w.r.t. time.

Model 1 - Implementation

Define mesh and spde model

border_simplified = st_simplify(border, dTolerance = 25)

mesh = fm_mesh_2d(boundary = border_simplified,
  max.edge = c(30, 90), cutoff = 15,
   offset = c(-0.05, -0.25),
  crs = st_crs(border))

spde = inla.spde2.pcmatern(mesh = mesh, 
  prior.range = c(200, 0.5), # P(range < 100) = 0.5
  prior.sigma = c(1, 0.5)) # P(sigma > 1) = 0.5

Note

  • Coordinate system is in Km (this is recommended always!)

  • The size of the domain is ca \(630\times488 \text{ Km}^2\) we use this as a guideline to define the prior for the range.

Model 1 - Implementation

Fit the model

cmp = ~ Intercept(1) + 
  time(time, model = "rw2") +
  space(geometry, model = spde)

lik = bru_obs(formula =  Y~ .,
              data = dd)

fit_1 = bru(cmp, lik)

# plot time effect

p = fit_1$summary.random$time %>% 
  ggplot() + 
  geom_ribbon(aes(ID, 
                  ymin = `0.025quant`,
                  ymax = `0.975quant`), 
              alpha = 0.5) +
  geom_line(aes(ID, mean))

Model 1 - Implementation

Space-time Predictions

pxl = fm_pixels(mesh, mask = border)

pxl_time = fm_cprod(pxl, data.frame(time = seq(1,12)))

pred_1 = predict(fit_1, pxl_time, ~ data.frame(mu = Intercept + space + time,
                                      time_effect = time,
                                      space_effect = space))

Spatio temporal models for areal data

The observation model

\[ Y_t(s)|\mu_t(s)\sim\mathcal{N}(\mu_t(s),\sigma_y^2); \qquad \eta_t(s) = \mu_t(s) \]

We are going to see 3 different models for the linear predictor \(\eta_{st}\)

  1. Model 1: \(\eta_t(s) = \beta_0 + u(s) + v_t\) with \(u(s)\sim\text{GF}(\rho,\sigma_u)\)

    Here the time effect \(v_t\) is constant w.r.t. space and the space effect \(u(s)\) is constant w.r.t. time.

  2. Model 2: \(\eta_t(s) = \beta_0 + u_t(s) + v_t\) with \(u_t(s)\sim\text{GF}(\rho,\sigma_u),\ t = 1,\dots,T\)

    Here spatial fields are different for each year, are independent and share range and sd.

Model 2 - Implementation

Fit the model

cmp = ~ Intercept(1) + 
  time(time, model = "rw2") +
  space(geometry, model = spde, 
        replicate = time)

lik = bru_obs(formula =  Y~ .,
              data = dd)

fit_2 = bru(cmp, lik)

# plot time effect

p = fit_2$summary.random$time %>% 
  ggplot() + 
  geom_ribbon(aes(ID, 
                  ymin = `0.025quant`,
                  ymax = `0.975quant`), 
              alpha = 0.5) +
  geom_line(aes(ID, mean))

Model 2 - Implementation

Space-time Predictions

pred_2 = predict(fit_2, pxl_time, ~ data.frame(mu = Intercept + space + time,
                                      time_effect = time,
                                      space_effect = space))

Spatio temporal models for areal data

The observation model

\[ Y_t(s)|\mu_t(s)\sim\mathcal{N}(\mu_t(s),\sigma_y^2); \qquad \eta_t(s) = \mu_t(s) \]

We are going to see 3 different models for the linear predictor \(\eta_{st}\)

  1. Model 1: \(\eta_t(s) = \beta_0 + u(s) + v_t\) with \(u(s)\sim\text{GF}(\rho,\sigma_u)\)

    Here the time effect \(v_t\) is constant w.r.t. space and the space effect \(u(s)\) is constant w.r.t. time.

  2. Model 2: \(\eta_t(s) = \beta_0 + u_t(s) + v_t\) with \(u_t(s)\sim\text{GF}(\rho,\sigma_u),\ t = 1,\dots,T\)

    Here spatial fields are different for each year, are independent and share range and sd.

  3. Model 3: \(\eta_t(s) = \beta_0 + u_t(s) + v_t\) with \[ u_t(s) = \phi\ u_{t-1}(s) + \omega(s),\text{ with } \omega(s)\sim\text{GF}(\rho,\sigma_u) \]

    Here we are adding more temporal dependence in the spatial term, using a AR1 type model Here the spatial fields are different for each year, they are independent on each other.

Model 3 - Implementation

Fit the model

h.spec <- list(rho = list(prior = 'pc.cor1', 
                          param = c(0, 0.9)))
cmp = ~ Intercept(1) + 
  time(time, model = "rw2") +
  space(geometry, model = spde, 
        group = time,
        control.group = list(model = 'ar1', 
                             hyper = h.spec))

lik = bru_obs(formula =  Y~ .,
              data = dd)

fit_3 = bru(cmp, lik)

# plot time effect

p = fit_3$summary.random$time %>% 
  ggplot() + 
  geom_ribbon(aes(ID, 
                  ymin = `0.025quant`,
                  ymax = `0.975quant`), 
              alpha = 0.5) +
  geom_line(aes(ID, mean))

Model 3 - Implementation

Space-time Predictions

pred_3 = predict(fit_3, pxl_time, ~ data.frame(mu = Intercept + space + time,
                                      time_effect = time,
                                      space_effect = space))

Adding covariates to space-time models

  • Covariates that only vary in space

    • Altitude, bathimetry,…

  • Covariates that vary in space and time

    • Temperature, precipitation, …

Space covariates

  • They should be stored in a terra, SpatRaster object.
library(terra)
library(tidyterra)
cov_space
class       : SpatRaster 
size        : 147, 147, 1  (nrow, ncol, nlyr)
resolution  : 4.465128, 3.11244  (x, y)
extent      : -465.4573, 190.9165, 7031.885, 7489.413  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=23 +south +datum=WGS84 +units=km +no_defs 
source(s)   : memory
name        :       last 
min value   : -3.1660518 
max value   :  0.8164945 
ggplot()  + geom_spatraster(data = cov_space) + geom_sf(data = dd) + scale_fill_scico()

Space covariates - Implementation

Setting up and running the model

cmp1 = ~ Intercept(1) +
  time(time, model = "rw2") +
  space(geometry, model = spde) +
  cov(cov_space,  model = "linear")

lik1 = bru_obs(formula =  Y~ .,
              data = dd_1)

fit_cov1 = bru(cmp1, lik1)

Look at the results

round(fit_cov1$summary.fixed,2)
           mean   sd 0.025quant 0.5quant 0.975quant  mode kld
Intercept  0.11 0.28      -0.45     0.12       0.66  0.12   0
cov       -1.42 0.09      -1.60    -1.42      -1.23 -1.42   0

Space-time covariates

Some covariates change in both space and time.

There are several ways of doing this….many quite confusing 🤪

Here is one recipe !

Space-time covariates

  • Have the covariate as a multilier spatraster object
# The covariate
cov_space_time
class       : SpatRaster 
size        : 147, 147, 12  (nrow, ncol, nlyr)
resolution  : 4.465128, 3.11244  (x, y)
extent      : -465.4573, 190.9165, 7031.885, 7489.413  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=23 +south +datum=WGS84 +units=km +no_defs 
source(s)   : memory
names       :          1,         2,         3,         4,         5,         6, ... 
min values  : -3.1660518, -1.275877, -1.572957, -1.991793, -2.704403, -3.526074, ... 
max values  :  0.8164945,  2.442684,  2.221178,  2.395833,  4.012937,  1.116201, ... 
# The data frame
dd_1[1:3,]
Simple feature collection with 3 features and 2 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -103.4191 ymin: 7460.922 xmax: -103.4191 ymax: 7460.922
Projected CRS: +proj=utm +zone=23 +south +datum=WGS84 +units=km +no_defs +type=crs
# A tibble: 3 × 3
      Y  time             geometry
  <dbl> <dbl>         <POINT [km]>
1  1.19     1 (-103.4191 7460.922)
2  2.72     5 (-103.4191 7460.922)
3  2.79     6 (-103.4191 7460.922)

Note The layers names in the raster are the same as the elements of the column time in the data.

Space-time covariates

Implementation

cmp1 = ~ Intercept(1) +
  time(time, model = "rw2") +
  space(geometry, model = spde) +
  cov(eval_spatial(cov_space_time, .data.,  time),  model = "linear")

lik1 = bru_obs(formula =  Y~ .,
              data = dd_1)

fit_cov1 = bru(cmp1, lik1)

round(fit_cov1$summary.fixed)
          mean sd 0.025quant 0.5quant 0.975quant mode kld
Intercept    0  0          0        0          1    0   0
cov         -1  0         -1       -1         -1   -1   0

The .data. indicates the data that are input in the bru_obs() function

Point Processes

Space-time Point Processes

  • For space-time point processes one can also use the group and replicate features

  • The one difference is that one has to define the integration scheme in both space and time!

Space-time Point Processes

  1. Define the component (with the group feature)
cmp <- ~ Intercept(1) +
  space_time(geometry,
             model = matern,
    group = season,
    ngroup = 4 )

lik =  bru_obs(
    geometry + season ~ .,
    family = "cp",
    data = mrsea$points,
    samplers = mrsea$samplers,
    domain = list(
      geometry = mrsea$mesh,
      season = seq_len(4)
    ))

fit <- bru(cmp, lik)

Space-time Point Processes

  1. Define the component (with the group feature)
  2. Define the likelihood model
cmp <- ~ Intercept(1) +
  space_time(geometry,
             model = matern,
    group = season,
    ngroup = 4 )

lik =  bru_obs(
    geometry + season ~ .,
    family = "cp",
    data = mrsea$points,
    samplers = mrsea$samplers,
    domain = list(
      geometry = mrsea$mesh,
      season = seq_len(4)
    ))

fit <- bru(cmp, lik)

Space-time Point Processes

  1. Define the component (with the group feature)
  2. Define the likelihood model
  • Where is our process defined (space and time)
cmp <- ~ Intercept(1) +
  space_time(geometry,
             model = matern,
    group = season,
    ngroup = 4 )

lik =  bru_obs(
    geometry + season ~ .,
    family = "cp",
    data = mrsea$points,
    samplers = mrsea$samplers,
    domain = list(
      geometry = mrsea$mesh,
      season = seq_len(4)
    ))

fit <- bru(cmp, lik)

Space-time Point Processes

  1. Define the component (with the group feature)
  2. Define the likelihood model

  • Where is our process defined (space and time)

  • Domain of integration

cmp <- ~ Intercept(1) +
  space_time(geometry,
             model = matern,
    group = season,
    ngroup = 4 )

lik =  bru_obs(
    geometry + season ~ .,
    family = "cp",
    data = mrsea$points,
    samplers = mrsea$samplers,
    domain = list(
      geometry = mrsea$mesh,
      season = seq_len(4)))

fit <- bru(cmp, lik)

Space-time Point Processes

  1. Define the component (with the group feature)
  2. Define the likelihood model
  3. Run the model (as usual)
cmp <- ~ Intercept(1) +
  space_time(geometry,
             model = matern,
    group = season,
    ngroup = 4 )

lik =  bru_obs(
    geometry + season ~ .,
    family = "cp",
    data = mrsea$points,
    samplers = mrsea$samplers,
    domain = list(
      geometry = mrsea$mesh,
      season = seq_len(4)))

fit <- bru(cmp, lik)

Space-time Point Processes

Results

ppxl <- fm_pixels(mrsea$mesh, mask = mrsea$boundary, format = "sf")
ppxl_all <- fm_cprod(ppxl, data.frame(season = seq_len(4)))

lambda1 <- predict( fit,ppxl_all,
  ~ data.frame(season = season, lambda = exp(space_time + Intercept)))

Non-separable models in inlabru

  • Some non-separable models are implemented in inlabru. See

Lindgrenn et al., A diffusion-based spatio-temporal extension of Gaussian Matern > fields (2024) SORT 48